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ABSTRACT 


A  simplified  technique  for  solving  "stiff"  differential  equations 
common  to  chemically  reacting  quasi- one -dimensional  inviscid  non¬ 
equilibrium  flow  analyses  is  presented.  The  method  is  basically  that 
used  by  Degroat  and  Abbett  with  the  generalization  made  to  allow  for 
a  prescribed  area  distribution  rather  than  the  constant  pressure  (and 
hence  constant  velocity)  process  assumed  by  them.  In  addition,  the 
equations  are  written  for  hydrogen/ air  combustion  rather  than  methane/ 
air.  The  advantage  of  the  technique  is  that  the  chemical  kinetic  calcu¬ 
lations  can  be  solved  quickly  without  loss  of  accuracy,  and  thus  can  be 
relegated  to  the  status  of  a  subroutine  in  complicated  fluid-dynamic 
problems. 
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SECTION  I 
INTRODUCTION 


Analysis  of  gaseous  flow  fields  in  which  the  gas  is  not  in  chemical 
equilibrium  consists  of  point  by  point  evaluation  of  the  gas  composition 
and  fluid- dynamic  variables.  Formal  analysis  involves  the  integration 
of  a  set  of  coupled,  first-order,  ordinary,  nonlinear,  stiff,  differential 
equations.  Numerical  solution  of  these  equations  often  eludes  the  com¬ 
puter  when  sought  by  classical  means  and  at  best  requires  considerable 
computation  time.  The  primary  reason  equations  of  this  class  (stiff) 
present  such  difficulties  is  that  an  exceedingly  small  integration  step- 
size  is  sometimes  required,  making  classical  integration  techniques 
impractical  even  on  the  most  sophisticated  computing  machines. 
Various  works  have  been  published  on  this  subject  (Refs.  1  through  7). 
The  method  of  solution  described  in  this  report  represents  a  combina¬ 
tion  of  some  of  the  improvements  made  in  several  of  the  above  refer¬ 
ences,  notably  Refs.  5  and  6.  In  this  method  the  stiffness  and  non¬ 
linearity  are  removed  entirely  by  using  approximation  techniques. 
Large  integration  stepsizes  become  possible  with  this  method,  particu¬ 
larly  for  near- equilibrium  flows. 


SECTION  II 
BASIC  EQUATIONS 


Let  y^  be  the  concentration  of  species  i  (moles/ cm^)  and  Mi  be  its 
molecular  weight.  Then 

yi  =  pT,  (1) 

where  ai  is  the  mass  fraction  of  species  i  and  p  is  the  mixture  density. 


At  any  instant  of  time,  the  rate  of  change  of  the  concentration  of 
species  i  in  reaction  j  is  given  by  the  general  formula: 


-4,7  « 

where  i/^and  are  either  zeros  or  integers,  and  the  factors  fj  and  bj 
are  forward  and  backward  reaction  rate  coefficients  for  the  jth  reaction, 
and  are  functions  of  temperature  only.  The  total  rate  of  change  of  the 
concentration  of  species  i  caused  by  all  the  chemical  reactions  is  given 
by  the  following  species  continuity  equation: 


yi  =  >'ij  =  I' i  CL  yi*  ya« "  “  “  ■ 
i=i 


(3) 


1 


AEDC-TR-68-268 


Suppose  that  in  an  integration  step,  the  variation  of  T  is  negligible, 
then  an  analytical  solution  to  Eqs.  (3)  can  be  obtained: 

yi  =  Yih)  (4) 


However,  Eqs.  (3)  are  nonlinear  and  thus  it  is  not  possible  to 
integrate  them  directly.  Now,  if  Eqs.  (3)  are  linearized  by  dropping 
higher  order  terms  which  may  be  assumed  to  be  negligibly  small  in  an 
integration  step,  then  an  exact  analytical  solution  can  be  obtained  over 
an  integration  interval.  This  has  been  done  (Ref.  5)  and  with  certain 
other  modifications  (described  in  detail  later)  the  concentration  can  be 
obtained  as  a  function  of  distance  (x)  also:  yj  =  yi(x). 

Since  the  solution  is  obtained  in  closed  form  over  an  integration 
interval,  no  stability  problems  arise.  The  advantage  of  this  technique 
stems  from  the  fact  that  the  integration  interval  stepsize  in  the  linear¬ 
ized  system  can  be  several  orders  of  magnitude  greater  than  in  the 
original  system.  In  the  next  sections,  a  specific  example  is  examined 
in  detail. 


SECTION  III 

COMBUSTION  OF  HYDROGEN  IN  AIR 


Let  air  be  represented  by  a  mixture  of  nitrogen  and  oxygen.  The 
nitrogen  is  considered  inert  and  the  oxygen  can  dissociate  or  react 
with  hydrogen.  Six  species  are  considered  in  the  chain  reaction  com¬ 
bustion  model:  H,  O,  H2O,  OH,  O2,  and  H2,  which  are  numbered  one 
through  six.  The  inert  nitrogen  is  numbered  seven  for  summation 
purposes  in  computing  the  third  body  M  for  third-order  reactions.  The 
reactions  considered  representative  for  combustion  of  hydrogen  in  air 
at  high  temperatures  are  (Ref.  5): 


Hj  +  0,  -  20H  j  =  1 

H  +  0a  fi  OH  +  0  2 

0  +  H,  «t  Oil  -  H  3 

Hj  +  OH  ^  H  +  1I30  4 

OH  +  OH  ^  0  +  H20  5 

H,  +  M  e  2H  +  M  6 

HjO  +  M  ^  OH  +  H  t  M  7 

OH  +  Mt»0  +  H  +  M  8 

02  +  M  ^  0  +  0  +  M  9 
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Equations  (2)  are: 


y4i  = 
y  22  = 
y  1 3  = 

>'l4  = 
y23  = 

Yu  = 

>17  = 

>'l6  = 

y»  = 


2fl)'s  >'6 

tzYiYs  -  b2y2y4 
fjyiVs  -  bjYjV* 
f4y4yt  -  b4ytyj 
Uy\  -  b5y2y, 
2fsysY  ■  2b6y  *Y 
f7ysY  -  b,y,y4Y 
f6y4Y  -  b,yiy2Y 
2f,y,Y  -  2b»y2Y 


(5) 


i  =  1 
H 


2 

0 


HaO  OH 


5 

02 


6 

H, 


i 

n2 


where 

v  -  E  »i 

i=i 

The  fj  and  bj  are  given  in  Appendix  II. 


Equations  (3) 

are  obtained  from  the  following: 

>n 

= 

0 

>21 

= 

0 

yn 

= 

0 

yu 

= 

-y22 

ya2 

= 

y>2 

>32 

= 

0 

yn 

= 

yn 

>23 

= 

-y» 

0 

>33 

• 

= 

0 

>14 

= 

>14 

>’24 

= 

ysx 

= 

yn 

v-x. 

= 

0 

>'« 

= 

Yi> 

yss 

= 

Yza 

yi« 

= 

yn 

>2  6 

= 

0 

yi* 

= 

0 

y  it 

= 

y  I? 

Yzi 

= 

0 

>'37 

= 

->'l7 

yu 

= 

>16 

Yi> 

= 

yi» 

yse 

= 

0 

ym 

= 

0 

• 

>29 

= 

>'29 

yss 

=■ 

0 

>'41 

= 

y4i 

>51 

= 

-y4./2 

hi 

= 

->'4l/2 

>42 

= 

>'22 

>52 

= 

~ >'22 

>62 

= 

0 

>43 

= 

>'l3 

>53 

= 

0 

>63 

= 

-yis 

>'44 

55 

->'14 

ys4 

= 

0 

ye4 

= 

-yi4 

y4s 

= 

-2y,s 

>'55 

= 

0 

yes 

= 

0 

>46 

=T 

0 

>'56 

= 

0 

y6« 

= 

-yn/2 

y47 

= 

y  17 

>'57 

= 

0 

>'67 

= 

0 

>'4. 

5 

->’i« 

>5. 

= 

0 

y66 

= 

0 

%* 

= 

0 

y*» 

= 

-y»/2 

ye9 

= 

0 

(6) 
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9 


1.6 


(7) 


The  following  two  equations  are  obtained  as  linear  combinations  of 
Eqs.  (7) 


2y«  +  yi  +  y4  +  2y,  =  o 


2yj  +  y 2  +  y*  +  y3  =  0 


(8) 


These  equations  express  the  conservation  of  the  number  of  atoms 
of  hydrogen  and  oxygen. 

The  following  additional  equations  must  be  added  to  make  the  problem 
determined. 


3.1  CONSERVATION  OF  SPECIES 


pa 


dctj 

77 


Mi  yi 


(9) 


where  d/  dx  means  substantial  derivative  and  u  is  the  flow  velocity. 


3.2  GLOBAL  CONTINUITY 


puA  =  m  =  constant 


(10) 


where  A  is  the  flow  area  and  m  is  the  mass  flow  rate. 

3.3  MOMENTUM  EQUATION 


dp 

dx 


+ 


P  U 


du 

dx 


0 


where  p  is  the  static  pressure. 


(ID 


3.4  ENERGY  EQUATION 


H„  =  II  + 


constant 


(12) 


4 


AEDC-TR-68-268 


where  H  is  the  static  enthalpy  of  the  gas  mixture  and  includes  sensible 
and  chemical  energies,  H0  is  the  total  enthalpy. 

H  =  £  as  hi  (13) 

1=1 

where  hj  is  the  sensible  enthalpy  of  species  i.  These  partial  enthalpies 
are  defined  through  partition  functions  of  the  species  which  are  functions 
of  temperature.  For  practical  purposes,  piecewise  parabolic  fits  can 
be  obtained  by  curve  fitting  data  such  as  are  tabulated  in  Ref.  8  thus: 

h;  =  A;  +  BjT  +  Cj  (T  -  Tn;)*  (14) 

where  Aj,  Bi,  and  Ci  are  polynomial  coefficients  and  T0^  is  a  reference 
temperature.  Details  are  given  in  Appendix  III. 


3.5  STATE  EQUATION 

7 

p  =  RT  £  y  i 
i=l 

where  R  is  the  universal  gas  constant. 


SECTION  IV 

LINEARIZED  EQUATIONS  FOR  THE  CONCENTRATIONS 


(15) 


The  mass  fractions  (<*i)  can  be  eliminated  by  means  of  Eq.  (1). 


**1 

yi  =  liT 

(16) 

yi  =  PYi 

(17) 

From  Eq. 

(9) 

<iy,  1 

—  =  —  ^ 
dx  p  u 

(18) 

From  Eq. 

(10) 

.  (I)''  y, 

dx  \A/ 

(19) 

In  every  computational  step,  the  following  simplifying  assumptions 
are  made: 

1.  The  reaction  rate  coefficients  fj  and  bj  are  constant 
throughout  a  step. 
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2.  The  logarithm  of  the  density  is  a  linear  function  of 
distance  throughout  a  step.  Consequently,  a  parameter 

8  =  -  *!.  (20) 
p  dx 

constant  throughout  a  step,  can  be  introduced.  Hence  if  p°  is 
the  initial  value  of  the  density  in  a  step  and  pc  is  the  final 
value,  the  ratio  ip  =  pc/p°  is  given  by 

ip  =  Pc./pc  =  1  +  <5Ax  (21) 

3.  The  sum  of  the  concentrations  Y  is  also  constant  through¬ 
out  a  step. 


Now  Eqs.  (5)  and  (7)  can  be  used  to  eliminate  the  terms  yi  in 
Eqs.  (19).  The  result  is  a  system  of  ordinary  differential  equations  in 
the  concentrations  which  can  be  solved  if  the  area  (A)  is  given  as  a 
function  of  x.  An  attempt  can  be  made  to  solve  these  equations  in  closed 
form  throughout  a  step. 


First,  Eqs.  (8)  can  be  used  to  eliminate  two  of  the  concentrations. 
The  slowly  varying  concentrations  Oo  and  H2  are  chosen.  Differentiating 
Eq.  (1)  yields 


dy*  Oi  ip  p  da* 
dx  Mj  dx  M{  dx 

and  using  Eq.  (9) 

dyt  1  dp  I  . 

dx  ~  J  ~ix  Yi  ^  u 

=  -  j'i  +  fyi 

u 


(22) 


(23) 

(24) 


Thus  using  Eqs.  (8)  and  (24)  one  can  obtain 


ys  =  - 


Yi 


y.  -  -  y‘  j  -  -  y,  +  *  jy?  +  y?  +  £ 


7*} 


where  ( )°  means  initial  value  at  the  beginning  of  the  step. 

L 


t  .  l  ,  /o  y°  +  y?  +  y?  I 
Let  b  =  p  !>p5  +  - - - | 

and  c  =  ip  +  y°  +  —  —  -  *  j- 


(25) 


(26) 
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Equations  {25)  and  (26)  can  be  restated: 

Yi  *  y3  +  >-4 


y  i  +  y« 

y«  =  c - - - y3 


(27) 


At  this  stage,  the  expressions  (Eqs.  (27))  for  y5  and  y0  can  be  sub¬ 
stituted  into  Eqs.  (5)  and  the  results  used  in  the  first  four  rows  of 
Eqs.  (7).  Therefore,  y^,  y2,  y3,  and  y^.,  now  depend  only  on  the 
variables  yj,  y2,  y3,  y4,  and  6  and  contain  a  number  of  constant 
coefficients  defined  by  the  values  of  T  and  of  each  concentration  at 
the  beginning  of  an  integration  step.  If  Eqs.  (19)  are  used  with  i  =  1, 

1  dp 

4;  6  =  —  — ,  the  yi  described  previously  and  a  known  area  distribution, 

a  system  of  four  ordinary  differential  equations  in  the  four  unknowns 
7i*  72*  73*  and  74  is  obtained  with  x  as  the  independent  variable. 

These  equations  have  the  general  form 


-5l  =  Fi(yk)  6, A)  i,k  =  1,  4  (28) 

dx 

where  the  Fi  (y&  6,  A)  are  linear  combinations  of  constant  terms,  first- 
order  terms,  and  second-order  terms  in  y^.  The  value  of  6  to  be  used, 
described  later,  does  not  complicate  the  system  since  it  is  assumed 
constant  throughout  a  step. 


A  system  such  as  Eqs.  (28)  cannot  be  solved  in  closed  form,  but 
can  be  easily  reduced  to  a  linear  system  by  writing 


y iy k  =  -  y i y k  +  yiyk  -r  ykyt 


(29) 


After  substitution  of  all  the  terms  yiyk  with  Eq.  (29),  the  system  of 
Eqs.  (28)  is  approximated  by  the  linear  system 


dx 


+ 


c; 


i  =  1,4 


(30) 


Appendix  IV  gives  the  expressions  for  the  constant  coefficients  a^k  and 


Several  methods  may  be  used  to  solve  Eqs.  (30). 
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4.1  METHOD  OF  MORRETTI  (REF.  5) 

In  this  analysis  of  hydrogen/air  combustion  the  velocity  was  not  a 
variable  and  hence  Eq,  (24)  is  solved  instead  of  Eq.  (30): 

— —  =  yi  +  5>'i  i  =  1,  4  (31) 

dt 

A  solution  to  this  set  of  equations  is 

Vi  -  Z  Ai.t  ="■*■  (32) 

k=l 

where  the  eigenvalues  r^  are  the  roots  of  the  characteristic  equation 

ai,k  -  rk  Si,k  =  0 

the  ^  are  the  corresponding  eigen  functions  and  6^  ^  is  Kronecker's 
delta.  In  general,  the  roots  may  be  complex  which  increases  the 
difficulty  of  obtaining  the  exact  solution  to  the  point  where  it  renders 
that  method  impractical  when  more  than  four  species  are  considered. 


4.2  METHOD  OF  DEGROAT  AND  ABBETT  (REF.6) 


Again  a  constant  velocity  system  was  considered  though  the  species 
considered  corresponded  to  combustion  of  methane.  An  approximate 
solution  to  Eq.  (31)  is  assumed 

yi  =  t  diq  t<i  (34) 

1=0 

This  solution  is  substituted  into  both  sides  of  Eq.  (31)  and  a  residue  is 
defined  as  follows: 

Rid)  =  -  n  -  Sy , 

=  £  qdi.it'1-  1  -  £  £  -  ci  (35) 

1=1  1=0  k=l 

Requiring  that 


yields  the  two  additional  conditions  required  for  the  evaluation  of  the 
diq,  d^q  for  each  equation 
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k=l 


n-  i 


it 

2 


n  A 

£  aik.dk  »t 
_k=  1 


0 


"A 


£ 

lr  =  l 


ajkdkqtg+l 
q  +  1 


Simplifying, 


aikdkot 


(37)* 


2 


I 


aikd|[Q 
q-  1 


0 


-  £  aikdko  (^j  -  c;  ^  -  0  (38) 

k=i 

where  is  the  number  of  rapidly  varying  species  (see  Eq.  (25)). 

There  results  a  set  of  2n^  linear  algebraic  equations  for  the  2n^ 
unknown  terms.  The  algebraic  equations  are  then  solved  by  a  maximum 
pivotal  point  matrix  reduction  method.  The  djq  being  known,  an  explicit 
solution  is  given  for  each  species  over  the  time  interval  considered. 


4.3  EXTENSION  OF  METHODS  1  AND  2  TO  VARIABLE  VELOCITY  SYSTEMS 


A  solution  to  Eq.  (30)  is  assumed  to  be  given  by 


2 


q=0 


Substitute  into  Eq.  (30) 


(39) 


(40) 


*The  inner  sumjnation  on  a^d^q  is  written  outside  the  brackets  in 
Ref.  6  which,  if  taken  literally,  would  lead  to  erroneous  results. 
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Differentiating  Eq.  (17) 

tyi  _  ^  rdyj 

dx  P  [_  dx 


A  new  residue  is  defined  as 


..  v  — 1  r  2  "A  q 

-  (t)  II  “ikdkqx'I  +  c5 
L9=°  k=l 


(41) 


(42) 


If  the  area  (A)  is  specified  as  a  function  of  x  in  the  interval,  then 
integration  of  the  residual  can  be  carried  out: 

Ax 

2  R,(x)  dx  =  0 


f 


Ax 

f 

J  Ax 


Ax 

2 


R;(x)  dx  =  0 


(43) 


and  again  a  set  of  2  n^  linear  algebraic  equations  is  obtained  which  can 
be  solved  by  standard  matrix  methods. 


The  unknown  parameter  5  is  required  to  enable  the  unknown  d^q 
terms  to  be  evaluated  via  Eqs.  (43).  There  is  no  way  in  which  6  can  be 
calculated  directly,  hence  an  iterative  solution  is  required.  A  simple 
test  is  used  to  verify  whether  the  assumed  value  of  6  is  compatible  with 
the  changes  in  species  concentration  calculated  from  a  solution  to 
Eqs.  (43).  The  temperature  at  the  end  of  the  integration  interval  can 
be  obtained  from  two  independent  equations:  (a)  equation  of  state,  and 
(b)  energy  equation. 


a.  The  pressure  (pc)  at  the  end  of  an  interval  is  obtained  from  a 
solution  of  the  continuity  and  momentum  equations 

jjp  _  2  /jn\2  P  dA  _  I  jjffl  (44) 

dx  p  \A )  La  dx  p  dxj 
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Since  A  is  specified  as  a  function  of  x  in  an  interval  Eq.  (23) 
can  be  integrated  after  making  the  substitutions 

p  =  p°  cSx  and  —  — ■  =  8 
p  dx 

From  the  equation  of  state 

Tc  =  Pc  (45) 


The  value  of  pc  is  obtained  from  integration  of  Eq.  (44). 
b.  The  enthalpy  at  the  end  of  an  interval  is  given  by 


Hc 


i=l 


From  Eq.  (14)  hj  =  Aj  +  BiT  +  C;(T  -  T0j)1 
From  the  continuity  equation 


m  =  puA  uc 


(46) 


(47) 


If  the  expressions  for  p  and  A  as  functions  of  x  are  inserted  in 
Eq.  (47)  and  the  value  of  uc  substituted  in  Eq.  (46)  a  quadratic 
expression  in  T  is  obtained: 


A*TCJ  -  2B*Td  +  c*  =  Ho  - 


(48) 


A*TCJ  -  2B“TC  +  C*  -  jH0  -  =  0  ■  (49) 

which  can  be  solved  for  Tc.,The  terms  A*,  B*,  and  C*  are  defined 
in  Appendix  III. 

The  values  of  Tc  obtained  in  Eqs.  (45)  and  (49)  should  agree.  If 
they  do  not  a  new  6  is  assumed  and  the  calculation  for  species 
composition  and  temperature  is  repeated.  A  Newton-Raphson 
iteration  procedure  for  compatible  values  of  6  and  T  has  been 
found  to  be  quite  adequate  and  usually  takes  no  more  than  three 
iterations  to  obtain  a  good  solution. 


Having  solved  for  the  new  concentrations  y\,  y2,  ys,  and  y4,  we 
obtain  the  new  values  of  ys  and  y6  from  Eqs.  (27). 
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4.4  COMPUTATION  OF  DENSITY 


Density  is  defined  by 

p  =  E  >;i  Mi  (50) 

i=l 


Now  the  mass  fraction  of  the  inert  species  (N2  =  07)  does  not  change; 
therefore. 


Pc  - 


(51) 


4.5  COMPUTATION  OF  VELOCITY 

A  new  velocity  is  obtained  from  the  continuity  equation 


using  the  correct  area  (Ac)  and  pc  from  Eq.  (51). 


4.6  COMPUTATION  OF  AREA 


The  flow  area  (A)  is  required  as  a  function  of  x  in  an  integration 
interval.  This  function  need  not  be  specified  as  a  function  of  distance 
from  the  origin;  thus,  any  complex  geometrical  configuration  can  be 
handled  readily  with  a  simple  function  for  the  integration  interval  only. 

Let 

Ac  =  A0  e£Ax  (52) 


where  Ac  is  the  area  at  the  end  of  an  interval,  and  AQ  is  the  area  at 
the  beginning  of  an  interval 

•••  e  -  E  (6  <53> 

Equations  (52)  and  (55)  are  to  be  used  in  Eqs.  (42),  (43),  (44),  and  (47) 

1.  If  the  area  is  constant  throughout:  j3  =  0. 

2.  If  the  pressure  is  constant  throughout:  j3  =  -6  (from 
the  continuity  equation). 

3.  If  the  pressure  is  specified  as  a  function  of  distance 
from  the  origin,  then  the  momentum  equation  and 
differential  form  of  the  energy  equation 


dH  do  _  _1  dp 

dx  dx  p  dx 


(54) 
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is  solved  for  Hc,  uc,  and  pc  and  the  area  can  be  obtained  from 
tfye  continuity  equation: 


A0 e0Ax  -  m 

(55) 

Pc  uc 

1  £»  r  m  i 

(56) 

Ao  Pc  ucJ 

4.7  STEPSIZE  CONTROL  AND  INITAL  STEPSIZE 


An  effective  stepsize  control  is  presented  based  on  the  error  gen¬ 
erated  in  the  solution  of  the  linearized  algebraic  equations.  This  is 
accomplished  by  a  comparison  of  the  results  obtained  from  the  two 
concentration  derivatives  (Eqs.  (40)  and  (41))  integrated  in  the  residual 
Eq.  (42). 

(yi)x  -  Ax  =  (yOx  +  Ay; 


Ay;  «  -H  .  AX 

dX 

Equation  40: 

(Ay;) i  =  i  (dj]  +  2d;2  AX  -  Sy;) 

i  =  l.nA 

Equation  41: 

(Ay;) 2  =  £  flikyk  -  c; 

'  V=1 

i  =  l>nA 

At  each  integration  step  a  comparison  is  made  of  the  above  calculations 
and  a  check  for  significance  of  the  error  is  also  included 


|(Ay;)2  -  (Ay;)j[ 
(yOx  +  Ax 


< 


significant  figure  in  the  value  of  (yk)x  +  AX 


SECTION  V 

NUMERICAL  CALCULATIONS  AND  RESULTS 


Method  3  has  been  programmed  and  found  to  be  a  very  effective 
method  of  analyzing  one -dimensional  (or  streamtube)  nonequilibrium 
flows.  The  calculated  results  are  at  least  as  accurate  as  those 
obtained  by  classical  methods,  and  computation  time  can  be  reduced 
drastically  in  the  cases  where  classical  methods  are  very  time 
consuming. 


5.1  EXAMPLE  OF  CONSTANT  PRESSURE  COMBUSTION 

This  example  involves  a  long  ignition  delay,  followed  by  a  short 
period  of  heat  release,  followed  by  along  period  of  near -equilibrium 
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flow.  A  Runge-Kutta,  predictor -cor rector  calculation  (Ref.  9)  was 
carried  out  for  a  duct  450  cm  long  having  an  initial  static  temperature 
of  1111°K  and  an  initial  static  pressure  of  0.  5  atm  in  2  hr  59  min  of 
computer  time  on  an  IBM  360/50.  The  results  are  shown  in  Pigs,  la 
through  e.  Appendix  I.  The  present  method  was  used  with  several 
values  of  stepsize  control  criteria:  0.0005;  0.005;  and  0.05.  No  sig¬ 
nificant  difference  occurred  between  the  Runge-Kutta,  predictor- 
corrector  results  and  the  present  method  except  for  the  largest  value 
of  stepsize  control  criteria  (0.  05)  where  differences  of  a  few  percent 
were  obtained  during  the  period  of  fast  heat  release.  Computation 
time  was: 

29  min  35  sec  at  0.0005  stepsize  control  criteria 
13  min  31  sec  at  0.  005  stepsize  control  criteria 
11  min  43  sec  at  0.  05  stepsize  control  criteria. 

The  savings  in  computation  time  were  appreciable.  The  differences 
of  a  few  percent  obtained  in  the  0.  05  control  case  are  somewhat  nebulous 
since  the  reaction  rate  data  are  not  known  sufficiently  accurately  to 
determine  quantitative  data  to  this  accuracy. 


5.2  EXAMPLE  OF  CONSTANT  AREA  COMBUSTION 

This  example  involves  a  short  ignition  delay  followed  by  a  period 
of  fast  reaction.  A  Runge-Kutta,  predictor- corrector  calculation  was 
carried  out  in  3  min  47  sec  for  a  duct  12.  7  cm  long  having  an  initial 
static  temperature  of  1160°K  and  an  initial  static  pressure  of  1.861  atm. 
The  results  are  shown  in  Figs.  2a  through  f.  The  present  method  was 
used  with  the  same  stepsize  control  criteria  with  no  significant  differ¬ 
ences  in  calculated  results.  Computation  time  was: 

3  min  48  sec  at  0.  005  stepsize  control  criteria 

1  min  59  sec  at  0.  05  stepsize  control  criteria. 

Again,  significant  savings  in  computation  time  were  obtained.  This 
particular  case  is  favorable  to  the  Runge-Kutta  predictor -corrector 
method  because  the  periods  in  which  the  equations  are  stiff  are  short 
(ignition  delay  and  approaching  equilibrium).  Further  savings  in  com¬ 
putation  time  are  obtained  if  certain  terms  in  the  linearized  coefficients 
(a^  k)  are  separated.  Only  a  small  fraction  of  the  a^  k  terms  involve 
the  parameters  6  and  0,  and  these  terms  only  need  be  calculated  at  each 
iteration. 
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5.3  EXAMPLE  OF  EXPANDING  STREAMTUBE  FLOW 
5.3.1  Specified  Area 


This  type  of  calculation  is  difficult  to  carry  out  for  two  reasons: 

1.  Species  production  equations  tend  to  "freeze"  as  the 
density  and  temperature  decrease. 

2.  The  approximate  Eq.  (30)  is  no  longer  a  linear  equation 

with  constant  coefficients  since  00  is  a  function  of 
streamtube  distance. 


Condition  (1)  tends  to  make  the  determinant  of  the  species  coef¬ 
ficient  matrix  tend  to  zero;  hence,  the  accuracy  in  theiinversion  of 
the  matrix  becomes  progressively  poor.  The  result  or  errors  caused 
by  condition  (1)  and  particularly  condition  (2)  appears  as  an  oscillation 
in  the  values  of  computed  variables.  The  amplitude  of  the  oscillation 
is  controlled  by  the  stepsize  control  criteria.  It  is  quite  small  for 
0.  005  as  shown  in  Figs.  3a  and  b.  The  geometry  is  shown  in  Fig.  3c, 
consisting  of  a  parabolic  turning  section  which  matches  the  slopes  of 
the  constant  area  combustor  and  a  15-deg  conical  expansion  section. 
The  oscillation  is,  of  course,  removed  by  setting  the  rate  Eqs.  (40) 
and  (41)  equal  to  zero  when  the  contribution  to  the  fluid-dynamic  equa¬ 
tions  is  reduced  to  a  negligible  amount  (flow  is  effectively  frozen). 


5.3.2  Specified  Pressure  as  a  Function  of  x 


The  integration  of  the  momentum  and  energy  equations  involves 
no  additional  errors  (the  area  equation  (Eq.  (52))  is  not  used)  and  no 
oscillation  is  obtained.  This  calculation  proceeds  much  faster  than 
the  above  specified  area  case. 

Since  the  specified  pressure  cases  generate  inherently  more  accu¬ 
rate  calculations,  it  is  preferable  to  evaluate  one  calculation  for  speci¬ 
fied  area  and  use  the  calculated  pressure  distribution  as  an  initially 
specified  function  of  x  for  additional  calculations.  This  is  also  suggested 
for  the  additional  reason  that  specified  area  distribution  cases  ignore 
wall  boundary  layers. 

5.4  FROZEN  FLOW 

This  calculation  required  about  20  sec  of  computing  time  for  the 
geometry  of  Fig.  3c.  The  species  rate  Eqs.  (40)  and  (41)  are  set  equal 
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to  zero  and  the  momentum  equation  is  integrated  directly  as  it  is  a 
function  of  the  assumed  density  change  and  given  area  variation  only. 
The  energy  equation  is  solved  as  before  to  obtain  a  temperature  which 
is  compared  with  the  value  from  the  equation  of  state.  An  upper  limit 
on  stepsize  is  necessary  to  ensure  smooth  plots.  A  temperature  and 
pressure  profile  is  shown  in  Figs.  4a  and  b  (data  correspond  to 
Section  III). 


5.5  EXAMPLE  OF  INFLUENCE  OF  INITIAL  CONCENTRATIONS 

The  Runge-Kutta  predictor-corrector  program  used  to  establish 
computer  time  for  comparison  purposes  employed  atomic  conservation 
equations  for  the  concentrations  of  hydrogen  and  oxygen  atoms.  Thus 
initial  values  of  these  concentrations  can  never  be  zero  since  round¬ 
off  error  in  the  initial  values  of  hydrogen  and  oxygen  molecules  gener¬ 
ates  initial  values  of  the  corresponding  atom  concentrations.  For  this 
reason  the  initial  values  generated  by  the  Runge-Kutta  predictor- 
corrector  program  were  used  as  initial  values  for  the  present  program 
and  are  shown  below  (mass  fractions  are  actual  input  values): 


i 

li 

moles /gm 

1 

(H) 

=  1.25992  x  10-8 

Density  (p)  =  5.  28314  x  10“4  gm/cm 

2 

(O) 

=  1. 96875  x 10'10 

Pressure  (p)  =  1.861  atm 

3 

(h2o) 

=  5. 55062  x  10"18 

Temperature  (T)  =  1159.  6°K 

4 

(OH) 

=  5.  87959  x  10"18 

Velocity  (u)  =  1.4129  x  108  cm/sec 

5 

(o2> 

=  6. 86381  x  10-3 

=  4635.  5  ft/sec 

6 

(h2) 

=  2.48016  x  10-3 

7 

<n2) 

=  2. 76755  x  IQ’2 

The  results  of  computations  with  these  values  were  discussed  earlier 
and  are  given  in  Figs.  2a  through  d. 

These  initial  values  for  hydrogen  and  oxygen  atom  concentrations 
completely  swamp  the  effect  of  including  reaction  1  as  suggested  in 
Ref.  10.  This  is  readily  illustrated  with  the  present  program  by  com¬ 
puting  the  results  with  these  two  concentrations  set  equal  to  some  negli¬ 
gible  number  (mass  fractions  of  lO-*8  were  used)  and  setting  reaction 
rate  f^  =  0.  These  results  are  shown  in  Figs.  5a  and  b.  If  the  initial 
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concentrations  are  maintained  at  this  negligible  quantity  and  reaction  rate 
fl  is  included,  the  results  are  again  similar  to  Figs.  2a  and  b  as  shown 
in  Figs.  5c  and  d.  Reaction  1  is  an  initiation  mechanism  which  is  neces¬ 
sary  to  describe  the  process  by  which  the  chain  reaction  ignition  reactions 
can  start. 

These  results  point  out  the  usefulness  of  being  able  to  set  the  initial 
concentrations  of  any  species  at  any  desired  value.  This  capability  is 
absolutely  essential  for  studies  of  ignition  delay  or  vitiation. 

SECTION  VI 

NORMAL  SHOCK  CALCULATIONS 


The  composition  of  the  streamtube  flow  is  assumed  to  be  frozen  just 
upstream  of  the  normal  shock.  The  values  of  fluid-dynamic  and  state 
variables  behind  the  normal  shock  are  obtained  by  an  iterative  solution 
to  the  conservation  equations  (Eqs.  (10),  and  (12)),  and  the  equation  of 
state  (Eq.  (15)).  These  computed  values  together  with  the  assumed 
frozen  composition  are  restored  to  the  nonequilibrium  program  as  initial 
values  and  the  approach  to  equilibrium  behind  the  normal  shock  is  then 
calculated.  A  normal  shock  at  any  station  in  the  combustor  or  nozzle 
can  be  evaluated. 
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APPENDIXES 

I.  ILLUSTRATIONS 

II.  REACTION  RATE  COEFFICIENTS 

III.  ENTHALPY  FIT  COEFFICIENTS 

IV.  LINEARIZED  COEFFICIENTS 
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a.  Concentration  Profiles 
Fig.  1  Constant  Pressure  Combustion 
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Concentration,  Log  yif  moles/cm3 
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b.  Concentration  Profiles 
Fig.  1  Continued 


Temperature, 
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b.  Density  Profile 
Fig.  2  Continued 
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c.  Nozzle  Geometry 
Fig.  3  Concluded 


-I _ I _ I _ I _ I  I  I 

32  36  40  44 


AEDC-TR-68-268 


Pressure,  atm 


Concentration,  Log  /] ,  moles/gm 
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Temperature, 
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Pressure,  atm 


5. 0  r  o  f  j  =  10^-  4  exp  (-  39, 000/RD  Atom  Mass  Fractions  Are: 
o  f !  -  1013- 4  exp  (-  39, 000/RD 

A  fj  »  0.0 

4.2  -  x  fj  -  0.0  For  this  Calculation, 

Initial  Atom  Mass 
Fractions  Are*. 

H-l.Z7xl0“8 

0-3.15xl0-9 


Log  (x  +  1) 


c.  Pressure  Profiles 
Fig.  5  Concluded 
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APPENDIX  II 

REACTION  RATE  COEFFICIENTS 


3  fj  Units 

1  1.0  x  1012-4  e"19-625/T  cm3/mole  sec 

2  3  x  10!4  e-8.  81/T 

3  3  x  10l4  e-4.  03/T 

4  3  x  10*4  e-3.  02/T 

5  3  x  10*4  e-3'  02/T 

6  1. 85  x  lOlfT-le-5*/?  cm6/mole2sec 

7  9.  66  x  lOiST-ie"62-  2/T 

8  8  x  1016  T_1  e-52.  2/T 

9  5.8  x  1016  T"1  e-60.6/T 

where  __ 

T  =  T/ 1000. 


bj  Ref. 

0  10 

2.  48  x  10^3e"®’  86/T  5 

1.3  x  10l4e-2.49/T 

1.  33  x  10l5e-10.  95/T 

3.  12  x  10^3e“^2.  51/T 
1016 

1017 
10*6 
6  x  10*4 
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APPENDIX  III 

ENTHALPY  FIT  COEFFICIENTS 


The  partial  enthalpies  (in  K  cal/gm)  are  fitted  by  means  of  the 
following  equations  (Ref.  5,  based  on  curve  fitted  data  from  Ref.  8): 


At  +  Bi T 

A;  +  B;T  +  Cj(T  -  Toi)J 

Di  -u  EiT 


■ 

Ai 

Bi 

Ci 

Di 

Ei 

T  • 
xoi 

1 

50.  22 

4.  93 

0.0 

0.0 

0.0 

6.0 

6.  0 

2 

3.  622 

0.3187 

0.0 

0.0 

0.0 

6.0 

6.0 

3 

-3.  3395 

0.4464 

0.0681 

-3.9456 

0.  7902 

0.5 

3.  94 

4 

0.4247 

0.4158 

0.0201 

0. 1631 

0. 5422 

0.5 

3.64 

5 

-0.0648 

0.2206 

0.0198 

-0. 2297 

0. 3168 

0.5 

2.  93 

6 

-1.004 

3.403 

0. 1968 

-4.286 

4.  831 

0.5 

4.096 

7 

-0.074 

0.  2488 

0.019 

-0. 1859 

0.3239 

0.  5 

2.48 

Coefficients  in  Eqs.  (48)  and  (49)  are: 


n’T 

A"  —  ^  Qj  Cj 

i=i 


tif 


1=1 


nf 

B*  =  £  «i(CiT0.  -  Bi/2)  C*  =  £  «i(Ai  +  Ci  Tqj) 

i=i 
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APPENDIX  IV 

LINEARIZED  COEFFICIENTS 


From  Eq.  (3) 


From  Eq.  (9) 


n  =  L  yjj 

j=l 


-©[I, alkyk  * C1] 

Using  Taylor  series  expansion  of  the  right-hand  side  of  Eq.  (IV-1): 

4  (5*°.. 

an  °Y  i) 


yij  -  y°i  +  l 


iyk  -  y°kf 


(IV-1) 

(IV- 2) 

(IV-3) 


k=i 

where  (  )°  means  evaluate  at  the  beginning  of  an  integration  interval. 
Equation  IV-1  becomes 

i'y°>  +  ih  fr lyk  -  yk' 


av-4) 


i=l 


k=i  j=i 


(IV- 5) 


and 


4 

-  £  aik  yk  +  Cj 

k=l 


(IV -6) 
(IV- 7) 
(IV- 8) 


(IV- 9) 
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The  yy  are  tabulated  as  follows: 


y» 

= 

0 

y»i 

= 

0 

Yu 

= 

0 

>41 

= 

>41 

y« 

= 

-Yu 

ha 

= 

yji 

y»2 

= 

0 

*42 

- 

y22 

y» 

= 

Yu 

Yu 

= 

-y.s 

Yu 

= 

0 

y42 

= 

Yu 

yu 

= 

y«4 

>'24 

= 

0 

yi4 

= 

yi4 

>44 

= 

— y  i4 

Yu 

= 

0 

Yu 

== 

yu 

Yu 

= 

Yu 

*45 

= 

-2y2> 

y«« 

= 

yu 

Yu 

= 

0 

Yu 

= 

0 

y* 

= 

0 

y» 

= 

Yu 

Yu 

= 

0 

ys? 

= 

-yn 

*47 

= 

yn 

yii 

= 

y» 

Yu 

= 

yis 

y»» 

0 

*45 

= 

-y» 

y»» 

= 

0 

Yu 

= 

y  29 

>'» 

= 

0 

= 

0 

(IV- 10) 


where  i  =  1,  2,  3,  and  4  are  the  only  species  considered  since  y$  and  yg 
are  given  by  Eq.  (25)  and  nitrogen  is  inert. 


Only  nine  of  the  terms  in  Eqs.  IV- 10  are  to  be  determined  since  the 
remaining  terms  are  either  zero  or  differ  by  some  integer  which  can  be 
positive  or  negative.  The  following  form  can  thus  be  used: 


where 


*ij 


4  9 

£  L  s»jRtj>rk  + 


k=lj=l 


4 


l 

If**' 


aik  =  L  SijRkj 
j=l 


C  =  l 
J= 


;  -  i  sijR°kj  ji 

l  L  k=l  J 


Si j  =  the  integers  in  Eqs.  (1V-10) 


Rkj 


yu 

3y°n 

dyk 


(IV- 11) 

(IV- 12) 

(IV- 13) 

(IV- 14) 
(IV- 15) 
(IV -16) 


The  R^j  are  determined  by  differentiating  Eqs.  (5)  after  substituion  of 
Eqs.  (25): 


y5  =  b  -  1/2  |ya  +  ys  +  y4l 
y«  =  c  -  1/2  iy,  +  2ys  +  y«i 


(IV- 17) 
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Substituting  in  Eqs.  (5) 

y«  =  2f,  [b  -  1/2  |y2  +  y,  +  y4fl  [c  -  1/2  {yx  +  2ys  +  y4  |] 

Yu  =  fayi  b  -  Y  Yi  hi  +  yj  +  )'«!  -  b2y,y4 

Yu  =  UYic - +  2y3  +  y4l  -  h,yjy4 

yi«  =  f«y«=  -  \  4  lyi  +  2ya  +  y4J  -  b4y,y2 
y«  =  fstf  -  bjyjyj 

yu  =  2f4Yc  -  f4Y  ly,  +  2ys  *  yJ  _  2b4y?  Y 
Yn  =  fryjY  -  b,yty4Y 


(IV- 18) 


y»  =  f.y-A'  -  b.y^jY 

y„  =  2f,Yb  -  f,Yly2  +  ys  +  y4!  -  2b,y’  Y 

Example: 

y  (yk  _  yok)  =  f2  [b  -  1/2  |y°  -  y°  +  y»n  (y,  -  y?)  -  P1  y?  +  bjy?l  (y2  -  y £ 
k=i  oyt  L2  j 

-  4  y°  (y>  -  yS)  -  [t  y?  +  b2  y°](y4  -  y?) 

=  Ri2(y,  —  y?)  +  (y*  —  y*)  +  Rj2(yj  -  y>)  +  ^«(y«  —  y«) 


The  coefficients  R^j  are: 


/ 

1 

2 

3 

1 

-ft  Con  * 

f2  Con  * 

- 

T  y^hsyll 

2 

-fj  Com  * 

-[r-  yi+b2^ 

J 

f3  Com  * 

3 

-fj  (Com+2Con)* 

T* 

- 

f3  y°2 

4 

-f^  (Com+Con)* 

“It  y^vi] 

- 

J  y°2+b3yi] 
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/ 

4 

5 

6 

1 

0 

-[f6Y+4W] 

2 

0 

-Vs 

0 

3 

-[f4y|+  b4y»] 

-  b5y°2 

-2f6Y 

4 

sfsyj 

"f6Y 

7 

8 

9 

1 

•v? 

-V°2Y 

0 

2 

0 

-b8yiY 

-[f9Y+4b9y“Y] 

3 

f7Y 

0 

-f9Y 

4 

-b7y°iY 

f8Y 

'f9Y 

where 

Con  =  b-1/2  jyjtyl+yjj 
Com  =  c-1/2  jyi+2y3+y4S 

By  inspection  only  seven  terms  involve  5,  j3  through  Con,  Com  and  only 
these  terms  (marked  *)  need  be  evaluated  at  each  iteration  for  6. 

The  following  coefficients  are  obtained  from  Eqs.  (5): 


R?i 

=  2hy*Y* 

Rf. 

=  2f6y°Y  - 

-  2b6y?3  Y 

R?, 

=  f3y?y° 

-  ivSW 

K 

=  f,y?Y  - 

bTyfy4°Y 

R?3 

=  - 

■  b,y?y? 

K 

=  f.y?Y  - 

b»y?y?Y 

r?< 

=  f4yJV?  - 

-  b4y?y? 

R?9 

=  2f,y?Y 

-  2b9y° J  Y 

r?5 

=  UyV  - 

•  bsy?y? 

43 


AEDC.TR-68.268 


44 


UNCLASSIFIED 

Security  Classification 


DOCUMENT  CONTROL  DATA  •  R  &  D 

(Security  classi  ficat  ion  of  title,  body  of  abstract  and  Indexing  annotation  must  be  entered  when  the  overall  report  Is  ctaaeltled) 


3  REPORT  TITL  E 

SIMPLIFIED  METHOD  FOR  SOLVING  PROBLEMS  INVOLVING  CHEMICALLY 
REACTING  ONE-DIMENSIONAL  FLOW 


4  OE5C  RIPTIVE  NOTES  (Type  ot  report  and  inclusive  dates) 

Final  Report  January  to  April  1967 

5-  AUTHOR(S)  (First  name,  middle  initial,  laat  name) 

I.  T.  Osgerby,  ARO,  Inc. 


10.  DISTRIBUTION  STATEMENT 

This  document  has  been  approved  for  public  release  and  sale; 
its  distribution  is  unlimited. 


13  ABSTRACT 

A  simplified  technique  for  solving  "stiff"  differential  equations 
common  to  chemically  reacting  quasi-one-dimensional  inviscid  non¬ 
equilibrium  flow  analyses  is  presented.  The  method  is  basically  that 
used  by  Degroat  and  Abbett  with  the  generalization  made  to  allow  for 
a  prescribed  area  distribution  rather  than  the  constant  pressure  (and 
hence  constant  velocity)  process  assumed  by  them.  In  addition,  the 
equations  are  written  for  hydrogen/air  combustion  rather  than  methane/ 
air.  The  advantage  of  the  technique  is  that  the  chemical  kinetic  calcu¬ 
lations  can  be  solved  quickly  without  loss  of  accuracy,  and  thus  can  be 
relegated  to  the  status  of  a  subroutine  in  complicated  fluid-dynamic 
problems . 


12.  SPONSORING  MILITARY  ACTIVITY 

Arnold  Engineering  Development 
Center  (AETS) ,  Arnold  Air  Force 
Station,  Tennessee  37389 


II  SUPPLEMENTARY  NOTES 

Available  in  DDC. 


fl  REPORT  DATE 

March  1969 

Bfl.  CONTRACT  OR  GRANT  NO 

F40600-69-C-0001 

b.  PROJECT  NO.  ' 

3012 

c-  Task  07 

d.  Program  Element  62402F 


7a.  TOTAL  NO  OF  PAGES  7b.  NO.  OF  REFS 

51  10 

9a.  ORIGINATOR'S  REPORT  NUMBERI5) 

AEDC-TR-68-268 

9b.  OTHER  REPORT  NO(S)  (Any  other  numbers  that  may  be  assigned 
this  report) 

N/A 


1  ORIGINATING  ACTIVITY  (Corporate  author) 

Arnold  Engineering  Development  Center 

ARO,  Inc. ,  Operating  Contractor 

Arnold  Air  Force  Station,  Tennessee  37388 


2«.  REPORT  SECURITY  CLASSIFICATION 


UNCLASSIFIED 


DD  ,fn°orv*\s1473 


UNCLASSIFIED 

Security  Classification 


Security  Classification 


